MODULE SRFWDIFSTL_MOD
CONTAINS
SUBROUTINE SRFWDIFSTL(KIDIA,KFDIA,KLON,KLEVS,&
 & PWSAM1M5 ,PCFW5  ,PRHSW5   ,PCDZ5 ,&
 & PWSADIF5 ,LDLAND ,LDALLAYS ,YDSOIL,&
 & PWSAM1M  ,PCFW   ,PRHSW    ,PCDZ  ,&
 & PWSADIF  )

USE PARKIND1  , ONLY : JPIM, JPRB
USE YOMHOOK   , ONLY : LHOOK, DR_HOOK, JPHOOK
USE YOS_SOIL  , ONLY : TSOIL

#ifdef DOC
! (C) Copyright 2012- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!**** *SRFWDIFSTL* -  DOES THE IMPLICIT CALCULATION FOR SOIL MOISTURE
!                     (Tangent linear)

!     PURPOSE.
!     --------
!          SOLVE TRIDIAGONAL SYSTEM OF EQUATIONS FOR SOIL MOISTURE.
!     IT SHOULD BE PRECEDED BY A CALL TO *SRFWEXC* AND FOLLLOWED BY
!     A CALL TO *SRFWINC*.

!**   INTERFACE.
!     ----------
!          *SRFWDIFSTL* IS CALLED FROM *SRFTSTL,SRFISTL*

!     PARAMETER   DESCRIPTION                                    UNITS
!     ---------   -----------                                    -----

!     INPUT PARAMETERS (INTEGER):
!     *KIDIA*      START POINT
!     *KFDIA*      END POINT
!     *KLON*       NUMBER OF GRID POINTS PER PACKET
!     *KLEVS*      NUMBER OF SURFACE LAYERS

!     INPUT PARAMETERS (LOGICAL):
!     *LDLAND*     LAND/SEA MASK (TRUE/FALSE)
!     *LDALLAYS*   TRUE - COMPUTATION FOR ALL LAYERS , FALSE - FOR TOP ONLY

!     INPUT PARAMETERS AT T-1 (REAL):
!  Trajectory  Perturbation  Description                                Unit
!  PWSAM1M5    PWSAM1M       PSI AT TIME LEVEL T
!  PCFW5       PCFW          MODIFIED DIFFUSIVITIES (LAMBDA-STAR IN DOCUM.)
!                            (INDEX JK REFERS TO LEVEL JK+1/2)
!  PRHSW5      PRHSW         RIGHT-HAND SIDE OF EQUATIONS
!                            TREATMENT OF EVAPORATION FLUX
!  PCDZ5       PCDZ          C*DZ  (PROFILE OF LAYER DEPTH FOR SOIL WATER)

!     OUTPUT PARAMETERS (REAL):
!  Trajectory  Perturbation  Description                                Unit
!  PWSADIF5    PWSADIF       PSI-STAR  DIVIDED BY ALFA

!     METHOD.
!     -------
!     *LU*-DECOMPOSITION AND BACK SUBSTITUTION IN ONE DOWNWARD SCAN
!     AND ONE UPWARD SCAN.

!     EXTERNALS.
!     ----------
!          NONE.

!     REFERENCE.
!     ----------
!          SEE SOIL PROCESSES' PART OF THE MODEL'S DOCUMENTATION FOR
!     DETAILS ABOUT THE MATHEMATICS OF THIS ROUTINE.

!     Original
!     --------
!     M. Janiskova              E.C.M.W.F.     08-02-2012

!     Modifications
!     -------------

!     ------------------------------------------------------------------
#endif

IMPLICIT NONE

! Declaration of arguments

INTEGER(KIND=JPIM), INTENT(IN)   :: KIDIA
INTEGER(KIND=JPIM), INTENT(IN)   :: KFDIA
INTEGER(KIND=JPIM), INTENT(IN)   :: KLON
INTEGER(KIND=JPIM), INTENT(IN)   :: KLEVS

LOGICAL,            INTENT(IN)   :: LDLAND(:)
LOGICAL,            INTENT(IN)   :: LDALLAYS

TYPE(TSOIL),        INTENT(IN)   :: YDSOIL

REAL(KIND=JPRB),    INTENT(IN)   :: PWSAM1M5(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PCFW5(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PRHSW5(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PCDZ5(:,:)
REAL(KIND=JPRB),    INTENT(INOUT):: PWSADIF5(:,:)

REAL(KIND=JPRB),    INTENT(IN)   :: PWSAM1M(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PCFW(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PRHSW(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PCDZ(:,:)
REAL(KIND=JPRB),    INTENT(INOUT):: PWSADIF(:,:)

!*         0.2    DECLARATION OF LOCAL VARIABLES.
!                 ----------- -- ----- ----------

REAL(KIND=JPRB) :: ZTCOE5(KLON), ZEBSW5(KLON,KLEVS)
REAL(KIND=JPRB) :: ZQDP5, ZDISC5,ZFAC5

REAL(KIND=JPRB) :: ZTCOE(KLON), ZEBSW(KLON,KLEVS)
REAL(KIND=JPRB) :: ZDISC, ZFAC, ZQDP, ZTPFAC2

INTEGER(KIND=JPIM) :: ILEVM1, JK, JL
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

!     ------------------------------------------------------------------

!*         1.    SET UP SOME CONSTANTS.
!                --- -- ---- ----------
IF (LHOOK) CALL DR_HOOK('SRFWDIFSTL_MOD:SRFWDIFSTL',0,ZHOOK_HANDLE)
ASSOCIATE(RSIMP=>YDSOIL%RSIMP)

ZTPFAC2=1.0_JPRB/RSIMP
ILEVM1=KLEVS-1

!*         1.1     SETTING OF RIGHT HAND SIDES.

DO JK=1,KLEVS
  DO JL=KIDIA,KFDIA
    IF (LDLAND(JL)) THEN
      PWSADIF (JL,JK) = ZTPFAC2*PWSAM1M (JL,JK)+PRHSW (JL,JK)
      PWSADIF5(JL,JK) = ZTPFAC2*PWSAM1M5(JL,JK)+PRHSW5(JL,JK)
    ENDIF
  ENDDO
ENDDO

!*         1.2     TOP LAYER ELIMINATION.

DO JL=KIDIA,KFDIA
  IF (LDLAND(JL)) THEN
    ZTCOE (JL) = PCFW (JL,1)
    ZTCOE5(JL) = PCFW5(JL,1)
    ZQDP5 = 1.0_JPRB/PCDZ5(JL,1)
    ZQDP  = -PCDZ(JL,1)*ZQDP5*ZQDP5
    ZDISC5 = 1.0_JPRB/(1.0_JPRB+PCFW5(JL,1)*ZQDP5)
    ZDISC  = -(PCFW5(JL,1)*ZQDP+ZQDP5*PCFW(JL,1))*ZDISC5*ZDISC5
    ZEBSW (JL,1) = (PCFW5(JL,1)*ZQDP5)*ZDISC+ZDISC5*PCFW5(JL,1)*ZQDP &
     & +ZDISC5*ZQDP5*PCFW(JL,1)
    ZEBSW5(JL,1) = ZDISC5*(PCFW5(JL,1)*ZQDP5)
    PWSADIF (JL,1) = ZDISC5*PWSADIF(JL,1)+PWSADIF5(JL,1)*ZDISC
    PWSADIF5(JL,1) = ZDISC5*PWSADIF5(JL,1)
  ENDIF
ENDDO


IF (LDALLAYS) THEN
!*         1.3     ELIMINATION FOR MIDDLE LAYERS.

  DO JK=2,ILEVM1
    DO JL=KIDIA,KFDIA
      IF (LDLAND(JL)) THEN
        ZQDP5 = 1.0_JPRB/PCDZ5(JL,JK)
        ZQDP  = -PCDZ(JL,JK)*ZQDP5*ZQDP5
        ZFAC  = ZTCOE5(JL)*ZQDP+ZQDP5*ZTCOE(JL)
        ZFAC5 = ZTCOE5(JL)*ZQDP5
        ZTCOE (JL) = PCFW (JL,JK)
        ZTCOE5(JL) = PCFW5(JL,JK)
        ZDISC5 = 1.0_JPRB/(1.0_JPRB+ZFAC5*(1.0_JPRB-ZEBSW5(JL,JK-1)) &
         & +PCFW5(JL,JK)*ZQDP5)
        ZDISC = -((1.0_JPRB-ZEBSW5(JL,JK-1))*ZFAC-ZFAC5*ZEBSW(JL,JK-1) &
         & +PCFW5(JL,JK)*ZQDP+ZQDP5*PCFW(JL,JK))*ZDISC5*ZDISC5
        ZEBSW (JL,JK) = (PCFW5(JL,JK)*ZQDP5)*ZDISC+ZDISC5*PCFW5(JL,JK)*ZQDP &
         & +ZDISC5*ZQDP5*PCFW(JL,JK)
        ZEBSW5(JL,JK) = ZDISC5*(PCFW5(JL,JK)*ZQDP5)
        PWSADIF (JL,JK) = (PWSADIF5(JL,JK)+ZFAC5*PWSADIF5(JL,JK-1))*ZDISC &
         & +ZDISC5*(PWSADIF(JL,JK)+ZFAC5*PWSADIF(JL,JK-1) &
         & +PWSADIF5(JL,JK-1)*ZFAC)
        PWSADIF5(JL,JK) = ZDISC5*(PWSADIF5(JL,JK)+ZFAC5*PWSADIF5(JL,JK-1))
      ENDIF
    ENDDO
  ENDDO

!*         1.4     BOTTOM LAYER ELIMINATION.

  DO JL=KIDIA,KFDIA
    IF (LDLAND(JL)) THEN
      ZQDP5 = 1.0_JPRB/PCDZ5(JL,KLEVS)
      ZQDP  = -PCDZ(JL,KLEVS)*ZQDP5*ZQDP5
      ZFAC  = ZTCOE5(JL)*ZQDP+ZQDP5*ZTCOE(JL)
      ZFAC5 = ZTCOE5(JL)*ZQDP5

      ZDISC5 = 1.0_JPRB/(1.0_JPRB+ZFAC5*(1.0_JPRB-ZEBSW5(JL,KLEVS-1)) &
       & +PCFW5(JL,KLEVS)*ZQDP5)
      ZDISC  = -((1.0_JPRB-ZEBSW5(JL,KLEVS-1))*ZFAC-ZFAC5*ZEBSW(JL,KLEVS-1) &
       & +PCFW5(JL,KLEVS)*ZQDP+ZQDP5*PCFW(JL,KLEVS))*ZDISC5*ZDISC5
      PWSADIF (JL,KLEVS) = (PWSADIF5(JL,KLEVS)+ZFAC5*PWSADIF5(JL,KLEVS-1)) &
       & *ZDISC &
       & +ZDISC5*(PWSADIF(JL,KLEVS)+ZFAC5*PWSADIF(JL,KLEVS-1) &
       & +PWSADIF5(JL,KLEVS-1)*ZFAC)
      PWSADIF5(JL,KLEVS) = ZDISC5 &
       & *(PWSADIF5(JL,KLEVS)+ZFAC5*PWSADIF5(JL,KLEVS-1))
    ENDIF
  ENDDO

!*         1.5     BACK-SUBSTITUTION.

  DO JK=ILEVM1,1,-1
    DO JL=KIDIA,KFDIA
      IF (LDLAND(JL)) THEN
        PWSADIF (JL,JK) = PWSADIF(JL,JK)+ZEBSW5(JL,JK)*PWSADIF(JL,JK+1) &
         & +PWSADIF5(JL,JK+1)*ZEBSW(JL,JK)
        PWSADIF5(JL,JK) = PWSADIF5(JL,JK)+ZEBSW5(JL,JK)*PWSADIF5(JL,JK+1)
      ENDIF
    ENDDO
  ENDDO
ENDIF


END ASSOCIATE
IF (LHOOK) CALL DR_HOOK('SRFWDIFSTL_MOD:SRFWDIFSTL',1,ZHOOK_HANDLE)

END SUBROUTINE SRFWDIFSTL
END MODULE SRFWDIFSTL_MOD
